AFRL-VS-TR-2001-1551 


Radiation  Belt  Test  Model 


John  W.  Freeman 


Rice  University 

Department  of  Physics  and  Astronomy 
6100  South  Main  St.,  Room  230 
Houston,  TX  77005 


October  2000 


Final  Report 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  IS  UNLIMITED. 


20020320  067 

AIR  FORCE  RESEARCH  LABORATORY 
Space  Vehicles  Directorate 
29  Randolph  Rd 

AIR  FORCE  MATERIEL  COMMAND 
Hanscom  AFB,  MA  01731-3010 


This  Technical  Report  has  been  reviewed  and  is  approved  for  publication. 


/I  yi^ 


CONTRACT  MANAGER 


GREGOKY/P.  GE^^ET,  Chief 
Space  Weather  Center  of  Excellence 


This  report  has  been  reviewed  by  the  ESC  Public  Affairs  Office  (PA)  and  is  releasable  to  the 
National  Technical  Information  Service. 

Qualified  requestors  may  obtain  additional  copies  from  the  Defense  Technical  Information 
Center  (DTIC).  All  others  should  apply  to  the  National  Technical  Information  Service  (NTIS). 


If  your  address  has  changed,  if  you  wish  to  be  removed  from  the  mailing  list,  of  if  the 
address  is  no  longer  employed  by  your  organization,  please  notify  AFRL/VSIM,  29  Randolph 
Rd.,  Hanscom  AFB,  MA  01731-3010.  ITtiis  will  assist  us  in  maintaining  a  current  mailing 
list. 

Do  not  return  copies  of  this  report  tmless  contractual  obligations  or  notices  on  a  specific 
document  require  that  it  be  returned. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


Public  reporting  burden  for  this  colieaion  of  information  is  esttmated  to  average  i  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewmg  the  colieaion  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspea  of  this 
collertion  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Direaorate  for  Information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202*4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduaion  Projert  (0704-0188),  Washington,  DC  20503. 


1.  AGENCY  USE  ONLY  (Leave  blank)  I  2.  REPORT  DATE 


Oct  2000 


4.  TITLE  AND  SUBTITLE 
Radiation  Belt  Test  Model 


6.  AUTHOR{S) 

John  W,  Freeman 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES} 

Rice  University 

Department  of  Space. Physics  and  Astronomy 
6100  South  Main  Street,  Room  230 
Houston,  TX  770C5 


3.  REPORT  TYPE  AND  DATES  COVERED 
FINAL  REPORT 


5.  FUNDING  NUMBERS 
PE:  62602F 
PR:  4026 
TA:  GC 
WU:  NA 


Contract:  F19628-97-C-001 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


9.  SPONSORING /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 
AFRL/VSBX 
29  Randolph  Road 
Ha-«scom  -B  MA  01731-3010 


10.  SPONSORING /MONITORING 
AGENCY  REPORT  NUMBER 


FRL-VS-TR-2001-1551 


C  ract  Mar  qer:  M.  HeinemL:  in/VSBX 


11.  SUPPLEMENTARY  NOUS 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 

Approved  for  Public  Release;  Distribution  Unlimited 


12b.  DISTRIBUTION  CODE 


13.  ABSTRACT  (Maximum  200  words) 

Rice  University  has  developed  a  dynamic  model  of  the  Earth's  radiation  belts  based 
.  on  real-time  data  driven  "boundary  conditions  and  full  adiabaticity.  The  Radiation 
'  Belt  Test  Model  (RBTM)  successfully  replicates  the  major  features  of  storm-time 
I  behavior  of  energetic  electrons:  sudden  -commencement  induced  main  phase  dropout  and 
recovery  phase  enhancement •  It  is  the  only  known  model  to  accomplish  the  latter. 

The  RBTM  shows  the  extent  to  which  new  energetic  electrons  introduced  to  the 
magnetosphere  near  the  geostationary  orbit  drift  inward  due  to  relaxation  of  the 
magnetic  field.  It  also  shows  the  effects  of  substorm  related  rapid  motion  of 
magnetotail  field  lines  for  which  the  3rd  adiabatic  invariant  is  violated.  The 
radial  extent  of  this  violation  is  seen  to  be  sharply  delineated  to  a  region  outside 
of  5Re,  although  this  distance  is  determined  by  the  Hilmer-Voigt  magnetic  field  model 
used  by  the  RBTM.  I 

The  RBTM  appears  to  provide  an  excellent  platform  on  which  to  build  parameterized  i 
refinements  to  compensate  for  unknown  acceleration  processes  inside  5Re  where 
adiabaticity  is  seen  to  hold.  Moreover,  built  within  the  framework  of  the  MSFM,  it 
offers  the  prospect  of  an  operational  forecast  model  for  MeV  electrons. 


14.  SUBJECT  TERMS  15.  NUMBER  OF  PAGES  | 

Radiation  belts,  magnetosphere,  energetic  electrons,  space  _ | 

weather  forecasting  16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 

UNCLAS 


18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 
UNCLAS 


19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 
UNCLAS 


20.  LIMITATION  OF  ABSTRACT 


NSN  7540-01-280-5500 


Standard  Form  298  (Rev.  2-89) 

Prescribed  by  ANSI  Std.  239-18 
298-102 


CONTENTS 


1.  Premise  and  objectives  of  the  contract 

2.  The  RBTM  1 

2. 1  Basics  of  the  model  1 

2.2  Installation  within  the  MSM  2 

2.3  Boundaries  and  support  parameters  2 

2.4  The  outer  boundary  flux  specification  4 

2.5  Preliminary  test  results  5 

3.  Delivered  product  9 

4.  Summary  and  recommendations  9 

4. 1  Summary  9 

4.2  Recommendations  10 

5.  Documentation  introduction  10 

Users  manual  and  code  description  1 1 

References  36 


iii 


1.  Premise  and  Objectives  of  the  Contract 

1.1  Premise 


The  premise  of  the  contract  was  that  all  adiabatic  behavior  of  energetic  electrons  must  be  examined 
in  a  fully  realistic,  time-variable  magnetic  field  before  the  effects  of  other  acceleration  or  radial  drift 
processes  can  be  fully  determined.  Once  adiabatic  effects  are  modeled  and  imderstood  the  amoimt 
of  additional  drift  or  acceleration  required  can  be  obtained  and  added  by  submodels  or  parameteriza¬ 
tion.  A  corollary  is  that  non-adiabatic  effects  can  be  detected  by  comparing  the  adiabatic  model  with 
satellite  observations.  Non-adiabatic  effects  may  be  treated  locally  (space  and  time-wise)  by  particle 
tracing  or  by  parameterization. 

1.2  Objectives 

Based  on  the  above  premise,  the  objectives  of  the  contract  were  to  1 .)  develop  an  energetic  electron 
radiation  belt  test  model,  RBTM,  based  on  the  fully  adiabatic  formalism  developed  by  Kim  and  Chan 
[1997];  2.)  test  the  model  against  satellite  data  and  then;  3.)  refine  the  model  for  improved 
performance  by  adding  new  drift  algorithms  or  parameterizations  that  would  correct  shortcomings 
turned  up  in  the  testing  process.  The  original  plan  was  for  phase  2,  the  testing,  to  be  performed  by  the 
Air  Force  Research  Laboratory.  The  last  two  phases  could  be  repeated  until  satisfactory  performance 
was  achieved.  Rice  would  develop  the  initial  model  and  perform  the  refinement  parameterization. 

2.  The  RBTM 


2. 1  Basics  of  the  model 


The  basic  engine  of  the  RBTM  is  the  fully  adiabatic  model  of  Kim  and  Chan  [1997].  This  model 
assumes  the  conservation  of  all  three  adiabatic  invariants,  which  amounts  to  the  assumption  that 
changes  in  the  magnetic  field  along  a  particle  path  occur  on  a  time  scale  that  is  long  compared  to  the 
characteristic  time  scale  of  each  of  the  three  invariants.  The  most  sensitive  of  these  is  the  third 
invariant  for  which  the  longitudinal  drift  time  for  a  1  MeV  electron  at  GEO  is  about  15  minutes.  The 
model  is  specialized  to  the  case  of  equatorially  mirroring  particles  so  that  the  second  invariant  is 
assumed  to  be  zero.  Using  Liouville’s  theorem  BCim  and  Chan  show  that  the  new  equatorial  flux  j, 
(electrons/unit  area,  time,  solid  angle  and  energy)  fix)m  time  1  to  time  2  is 


(1) 


where  the  E  is  the  kinetic  energy,  L  is  the  drift  shell  and  B}(JLj)  and  B2{L2)  are  the  magnetic  field 
strengths  on  the  respective  drift  shells  before  and  after  the  step. 


The  new  drift  shells  are  located  numerically  by  integrating  the  polar  flux  around  a  contour  of 
constant  B  imtil  a  drift  path  that  conserves  the  flux  enclosed  by  that  L  drift  shell  has  been  found. 
That  sets  the  new  B  for  the  determination  of  the  flux  according  to  equation  (1). 

The  new  energy  E2  is  determined  using  conservation  of  the  (relativistically  correct)  first  adiabatic 
invariant: 
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El  =  -fnc  +  ^E^  -  2mcE, )  +  (mc^f 


(2) 


The  actual  computational  logic  flow  of  the  RBTM  employs  the  following  steps: 

1.  Setup:  For  a  given  drift  shell  L„,  calculate  the  enclosed  ionospheric  polar  flux  using  Roederer’s 
expression  =  2ns^Bo/L„,  where  a  is  the  radius  of  the  Earth,  Roederer  [1970]; 

2.  Find  a  starting  equatorial  Btest  for  the  same  L„  using  the  dipole  expression  Btes^Bo/Ln  ; 

3.  For  each  local  time  grid  point  (j)  find  the  colatitude  for  using  the  Hilmer-Voigt  model; 

4.  Find  the  polar  flux  O/est  that  is  enclosed  by  these  colatitude  points  using  <t>,csr  =  2j  In  a^Bo  sin^ 
(colat  (i,j))  /  48; 

5.  Compare  On  and  0,ej,.  If  absolute  difference  >  0.01  then  iterate  using  a  new  Btest  until  satisfied; 
(The  new  B  is  foimd  by  a  bounded  search.) 

6.  Compute  new  flux  and  energies  at  each  grid  point  using  the  Kim  and  Chan  flux  and  energy 
advancement  equations  shown  above; 

7.  Average  the  computed  fluxes  over  the  drift  shell  and  interpolate  to  the  grid  points. 

2.2  Installation  within  the  MSM 

Early  in  the  contract  it  was  decided  to  install  the  RBTM  as  an  new  module  within  the 
Magnetospheric  Specification  Model,  MSM/MSFM.  This  approach  offered  several  advantages.  It 
would  provide  the  RBTM  access  to  the  Hilmer-Voigt  magnetic  field  model  within  the  MSM.  It 
would  provide  access  to  all  necessary  input  parameters.  It  would  allow  the  use  of  the  MSM  graphic 
output  for  the  energetic  electrons,  and  finally  it  would  permit  easy  transition  to  operation  since  the 
computer  interface  of  the  RBTM  would  be  imchanged  fi’om  that  of  the  MSM.  In  fact,  the  original 
MSM  lower  energy  electron  and  ion  species  output  would  still  be  available.  The  range  of 
applicability  of  the  MSM  would  be  extended  to  5  MeV  electrons. 

2.3  Boimdaries  and  support  parameters 

The  fully  adiabatic  model  is  a  boxmded,  gridded  model  that  requires  initial  electron  fluxes  at  all  grid 
points  at  the  time  of  a  cold  start  and  dynamic  input  fluxes  at  the  inner  and  outer  boundaries.  It  was 
clear  that  the  CRRESELE  statistical  model  could  be  used  to  provide  the  starting  electron  fluxes. 
CRRESELE  could  also  be  used  for  the  dynamic  fluxes  on  the  iimer  boundary  at  2.5  Re  provided  the 
time  resolution  of  CRRESELE,  one  day,  could  be  improved  to  match  the  cadence  of  MSM  with  its 
15  minute  time-step. 

The  outer  boundary  presented  a  more  difficult  challenge  because  of  known  local-time  variations  in 
the  energetic  electron  flux.  CRRESELE  is  symmetric  in  local  time  and  could  not  provide  the 
required  local-time  dependent  fluxes.  Local-time  variations  were  deemed  to  be  less  important  at  the 
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The  outer  boundary  presented  a  more  difficult  challenge  because  of  known  local-time  variations  in 
the  energetic  electron  flux.  CRRESELE  is  symmetric  in  local  time  and  could  not  provide  the 
required  local-time  dependent  fluxes.  Local-time  variations  were  deemed  to  be  less  important  at  the 
inner  boundary.  However,  Rice  had  experience  with  artificial  neural  network  systems  that  could  be 
trained  to  provide  the  continuous  dynamic  fluxes  on  a  model  outer  boundary  set  at  the  geostationary 
orbit.  The  neural  network  could  be  driven  by  solar  wind  data  only.  Moreover,  if  those  data  were 
obtained  from  an  LI  spacecraft  the  Model  would  have  intrinsic  forecast  capability. 

The  modeling  region  is  the  magnetic  equatorial  plane.  Figure  1  shows  the  basic  concept  of  the 
RBTM  with  the  modeling  region  and  the  initial  and  boundary  conditions  specified. 


The  RadiationBelt  Test  Model 


.  CRRESELE  provides  the 

Geostationary  Orbit  Initialization  flux  at  each 

Boundary  flux  specified  by  Solar  Wind  Grid  point 

driven  neural  networks 


Figure  1 .  Modeling  region  and  initial  and  boundary  conditions  for  the  RBTM 

Figure  2  is  a  flow  chart  showing  the  modules  of  the  revised  MS(F)M  as  envisioned  for  the  final 
product.  Solar  wind  data,  preferably  from  an  LI  spacecraft  is  fed  into  the  data  conditioning 
routines  within  the  MSM  and  into  KIPNET.  KPNET  is  used  to  generate  a  pseudo  Kp  at  15- 
minute  intervals.  This  is  converted  to  Apl5  which  is  then  used  to  generate  a  pseudo  Ap  15  at  15- 
minute  intervals  in  time  synchronism  with  the  MSM  time  steps.  CRRESELE  can  then  provide 
the  starting  and  inner  boundary  fluxes  to  the  Fully  Adiabatic  Model  (FAM).  Meanwhile  MSM 
is  providing  B-field  values  to  FAM  at  the  same  cadence  and  the  neural  network  using  data  from 
MSM  can  provided  the  outer  boundary  electron  fluxes.  (Note:  in  the  delivered  version  this 
function  is  replaced  by  an  interpolation  of  GEO  satellite  data.)  The  FAM  computed  fluxes  are 
available  for  output  in  the  MSM  graphic  display  system. 

An  operator  switch  (not  shown  in  Figure  2)  permits  the  electron  flux  output  to  default  to  the  pure 
CRRESELE  output,  thus  providing  a  15-minute  CRRESELE  update. 
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Radiation  Belt  Test  Model 

Flow  Chart  rbtm 

Output 
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Version  6 
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Figure  2.  Flow  chart  for  the  RBTM. 

2.4  The  outer  boundary  flux  specification 

The  proposed  plan  was  to  use  an  improved  version  of  a  neural  network  specification  algorithm 
produced  under  a  separate  contract  by  O’Brien  and  Freeman  [1998]  to  specify  energetic  electron 
fluxes  on  the  outer  boundary  (6.6  Re).  Upon  examination  it  was  determined  that  several  errors  that 
could  affect  the  accuracy  had  been  made  in  the  preparation  of  that  algorithm  and  that  the  neurons 
should  be  retrained.  Work  was  begun  on  this  diiring  the  summer  of  1999  by  an  imdergraduate  intern 
Cadence  Ellington.  Most  of  the  summer  effort  involved  setting  up  the  training  data  sets,  which  had  to 
be  rebuilt  jfrom  scratch.  As  a  result  the  retraining  was  not  completed  during  the  summer.  It  had  been 
intended  to  finish  that  job  during  the  summer  of  2000.  In  the  meantime,  coding  was  completed  by 
Bonnie  Hausman  on  the  rest  of  the  model  including  the  nemal  network  to  compute  the  15 -minute 
Apl5  to  drive  CRRESELE.  Since  the  GEO  electron  neural  network  was  not  ready,  and  an  outer 
boundary  energetic  electron  specification  was  needed,  it  was  decided  to  proceed  in  two  intermediate 
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steps.  First  we  used  CRRESELE  as  the  outer  boundary  specification.  This  did  not  provide  local  time 
variation  in  the  electron  flux  but  it  did  allow  us  to  begin  an  initial  debugging  of  the  code. 

Next  we  prepared  a  GEO  orbit  electron  flux  specification  for  the  November  3-8,  1993  magnetospheric 
storm  using  a  two-dimensional  interpolation  of  data  from  the  LANL  satellite  1990-046.  The  inter¬ 
polation  spread  the  data  in  local  time  and  imiversal  time  using  the  MATLAB  function  GRIDDATA. 
This  provided  a  useful  dynamic  test  input  data  set  for  the  GEO  electron  flux  with  the  required  local 
time,  including  diurnal,  variation.  This  flux  distribution  is  shown  in  Figure  3. 


1989-046  1  MeV  log  flux  w/cubic  interpolation  to  all  local  times 


Figure  3.  Plot  of  the  1  MeV  electron  data  from  satellite  1989-046  interpolated  in  Universal  Time 
and  Local  Time  as  used  as  input  for  the  outer  boundary  dynamic  flux  specification  in  the  RBTM. 

2.5  Preliminary  Test  Results 

Because  it  is  built  as  a  module  within  the  MSFM,  the  RBTM  uses  the  MSM  output  graphics  as  the 
primary  display  mode.  The  only  change  from  the  standard  MSM  display  was  to  cut  off  the  outer 
boundary  of  the  plotted  energetic  electron  flux  2  Re  beyond  GEO.  The  displayed  fluxes  between  GEO 
and  the  display  boundary  are  obtained  by  an  extrapolation. 

Figure  4  is  an  example  of  the  output  plot  from  the  November  1993  storm  run.  The  circle  just  inside 
the  model  outer  boundary  represents  the  geostationary  orbit.  The  Dst  curve  for  the  storm  is  also 
shown. 
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Exhaustive  testing  of  the  RBTM  against  satellite  data  was  intended  to  be  done  by  APRL  and  is 
therefore  outside  the  scope  of  the  contract,  however,  we  have  conducted  some  limited  tests  using  the 
November  1993  storm. 

Figure  5  shows  a  sequence  of  frames  throughout  this  storm  for  1.6  MeV  electrons  as  modeled  by  the 
RBTM.  The  frame  letters  refer  to  the  times  shown  in  the  Dst  curve  in  Figure  4.  Following  two 
prestorm  frames,  A  and  B,  the  compression  from  the  sudden  commencement  can  be  seen  in  frame  C. 
Frames  D  and  E  show  the  typical  depletion  at  GEO  associated  with  the  main  phase  of  the  storm.  The 
effects  of  a  substorm  can  be  seen  in  frames  D  and  E.  Dipolarization  of  the  tail  field  reduces  the  flux  at 
GEO  substantially.  This  flux  returns  as  the  tail  field  expands  outward  again  as  seen  in  frame  F. 
Frames  G  and  H  show  the  typical  recovery  phase  enhancement.  Note  that  the  outer  region  fluxes 
exceed  prestorm  values. 

In  these  runs  the  electron  fluxes  have  been  averaged  around  a  drift-shell  at  the  end  of  each  time  step 
thus  washing  out  all  local  time  variations  except  those  due  to  the  asymmetry  of  the  magnetic  field. 
This  was  done  to  maintain  consistency  with  the  full  adiabatic  concept  of  the  model.  However,  RBTM 
test  runs  not  incorporating  the  final  drift-shell  average  show  significant  local  time  variations  due  to 
rapid  field  line  motion  in  the  tail.  This  field  line  motion  occurs  between  15-minute  time  steps  of  the 
model  and  on  a  time  scale  of  the  order  of  the  electron  longitudinal  drift  time  thus  resulting  in  a 
violation  of  the  3’’*^  invariant  which  shows  up  in  these  (instmctive)  non-average  runs.  The  delivered 
version  of  the  RBTM  has  the  final  drift-shell  average. 

To  summarize,  these  preliminary  test  runs  indicate  that  the  model  is  functioning  as  intended  at  this 
stage. 

Based  on  the  foregoing  and  additional  inspection  of  the  model  output  compared  with  the  magnetic 
field  model,  we  can  make  the  following  statements  regarding  the  function  of  the  RBTM: 

1.  The  model  replicates  the  main  phase  dropout. 

2.  There  are  strong  perturbations  of  the  magnetic  field  in  the  midmght  region  just  outside  of  about  5 
Re  during  the  main  phase/substorm  period  that  break  the  3^*^  adiabatic  invariant  of  the  MeV  electrons. 

3.  The  model  reacts  to  these  perturbations  appropriately  by  pushing  the  electrons  inward  upon 
compression  and  pulling  them  outward  during  expansions,  however,  adiabaticity  caiuiot  be  trusted  at 
these  times.  A  more  elaborate  particle  tracing  scheme  is  needed. 

4.  There  is  a  relatively  sharp  inner  boundary  to  these  excursions  of  the  field  that  lies  just  inside  the 
geostationary  orbit;  fluxes  are  far  more  stable  inside  this  boundary. 

5.  Inside  5Re,  throughout  the  storm,  adiabaticity  probably  works  as  a  baseline  but  parameterization 
may  be  required  for  supplemental  acceleration. 

6.  During  the  recovery  phase  electrons  injected  at  GEO  propagate  inward  to  fill  the  region  into  at  least 
5  Re.  This  provides  some  recovery  phase  enhancement  in  the  region. 

7.  A  careful  comparison  with  satellite  data  is  needed. 

8.  A  phase  space  density  analysis  would  be  useful. 
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The  November,  1993  Storm  Dst 


Figure  4.  A  single  frame  example  of  an  output  plot  for  the  log  flux  of  1.6  MeV  electrons  for  09:15 
UT,  November  4,  1993.  The  Sun  is  on  the  left  and  the  tic  marks  are  5  Re.  The  circle  is  the 
geostationary  orbit.  The  color  bar  is  the  log  flux  in  electrons/cm^  s  sr  keV.  This  frame  occurs  just 
after  a  substantial  expansion  of  the  tail  field  following  a  substorm.  The  Dst  plot  is  also  shown. 
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Figure  5.  Sample  output  plots  of  the  RBTM  with  the  data  shown  in  Figure  3.  used  as 
input  for  the  outer  boundary. 
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3.  Delivered  product 

The  delivered  model  is  configured  to  run  operationally  using  geostationary  satellite  data.  These 
data  are  interpolated  in  local  time  and  universal  time  to  provide  the  real-time,  dynamic  outer 
boimdary  conditions  to  run  the  model.  This  mode  is  imtested  due  to  lack  of  time. 

The  model  will  have  two  output  modes: 

•  Standard  mode  in  which  the  output  will  be  the  RBTM  computed  electron  fluxes 

*  CRRESELE  output  in  which  tiie  fast  CRRESELE  fluxes  replace  the  standard  output.  This 
mode  is  operator-selectable  by  coding  the  ENCHAN  file  to  zero  energy  channels. 

The  delivered  model  uses  the  MSFM  interface  and  requires  only  solar  wind  and  LANL  data  to 
operate. 

4.  Summary  and  Recommendations 
4.1  Summary 

The  RBTM  runs  and  provides  a  useful  dynamic  simulation  of  energetic  electrons  in  the 
equatorial  plane  of  the  Earth’s  magnetosphere.  The  RBTM  has  not  been  tested  against  satellite 
data  and  certainly  needs  additional  algorithms  to  fine-tune  the  flux  values  and  energies  for 
greater  fidelity.  The  limited  testing  that  has  been  done  using  the  November  1993  storm  suggests 
the  following  observations: 

1 .  Rapid  storm/substorm  field  line  motion  beyond  about  5  Re  breaks  the  3^*^  adiabatic  invariant 
during  certain  times.  Development  of  supplemental  algorithms  to  address  this  region  at  these 
times  is  a  next  logical  step  in  the  refinement  of  this  code  but  is  beyond  the  scope  of  this  contract. 
Full  particle  tracing  may  be  required  for  this  region  and  for  certain  times,  with  a  “smart”  sensing 
algorithm  to  determine  when  the  full  trace  should  kick  in. 

2.  Inside  5  Re  the  fully  adiabatic  model  could  form  a  baseline  algorithm  upon  which  to  build 
with  parameterized  adjustments  to  the  flux  and  energy.  The  Kim-Chan  equations  (1)  and  (2) 
provide  a  convenient  fi'amework  for  accomplishing  this.  Coefficients  could  be  added  to  the  ratio 
of  B  terms  that  represent  the  rate-of-change  of  the  field.  Again,  this  cannot  be  done  imder  the 
present  contract. 

3.  Without  extensive  comparison  with  satellite  data  inside  geostationary  orbit,  it  is  not  possible  to 
determine  how  much  radial  motion  of  the  drift  shells  is  being  provided  by  the  storm-time  motion 
of  field  lines  and  therefore  how  much  radial  diffusion  would  need  to  be  added  to  the  model. 
Preliminary  inspections  suggest  that  the  region  just  inside  geostationary  orbit  fills  in  quickly 
during  the  recovery  phase  when  flux  enhancements  occur  at  GEO.  However,  this  might  be  an 
artifact  of  the  current  version  of  the  model  which  cannot  be  considered  to  be  fully  debugged. 
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4.  The  delivered  version  of  the  code  takes  averages  of  the  flux  around  each  newly  computed  drift 
shell.  This  is  done  to  preserve  the  concept  of  full  adiabaticity,  however,  test  versions  of  the  code 
that  did  not  that  take  T-shell  averages  revealed  interesting  local-time  asymmetries  in  the 
midnight  region  that  reveal  the  violation  of  the  flux  invariant. 

5.  The  concept  of  using  dynamic  satellite  data  at  GEO  as  an  outer  boundary  specification  along 
with  CRRESELE  for  the  inner  and  initial  specification  appears  to  work  very  well.  Adding  a 
neural  network  specification  for  the  outer  boundary  flux  would  free  the  model  from  1  AU  input 
data  and  convert  the  RBTM  to  a  true  forecast  model. 

6.  In  summary,  the  RBTM  concepts  have  been  shown  to  work  and  the  approach  has  been 
justified. 

4.2  Recommendations 

We  recommend  that  this  model  be  brought  to  completion  by  implementation  of  the  following 
steps: 

1.  Train  neural  networks  to  provide  the  geostationary  outer  boundary  fluxes. 

2.  Test  the  model  thoroughly  against  satellite  data  and  determine  areas  where  adiabaticity  is 
broken  and  acceleration,  loss  processes  and/or  radial  diffusion  is/are  needed. 

3.  Add  parameterization  or  individual  particle  drift  algorithms  to  overcome  the  shortcomings 
found. 

4.  Retest  and  adjust  parameterizations  for  optimal  accuracy. 

5.  Replace  existing  MSM  with  this  upgrade  in  operations. 

5.  Documentation 

The  following  pages  describe  the  RBTM  code.  The  MSFM  is  not  discussed  here  since  there  exists  a 
separate  final  report  for  MSFM  delivered  to  the  Air  Force  on  February  26,  1994  under  contract 
F19628-90-K-0012.  The  output  and  input  formats  have  not  changed.  All  previous  programs  and 
graphics  routines  developed  for  the  MSFM  will  work  with  the  RBTM. 
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RBTM  Users  Manual 


Types  of  input  files  for  RBTM: 

1)  static  files  required  for  run 

2)  data  files  required  for  run 

3)  data  files  optional  for  run 
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Standard  input  format  for  RBTM  subroutine  INDATA 

Example 

4.000  1993.000  2.333  -999  -999  -999  -999  -999 

4.333  1993.000  2.666  -999  -999  -999  -999  -999 

5.000  1993.000  2.999  -999  -999  -999  -999  -999 

The  8-word  record  is  as  follows: 

Word  1  =  data  value  at  t 
Word  2  =  year  at  t  (4-digits) 

Word  3  =  decimal  day  at  t 

Word  4  =  spacecraft  geomagnetic  latitude  at  t 

Word  5  =  spacecraft  geomagnetic  longitude  at  t 

Word  6  =  spacecraft  altitude  at  t 

Word  7  =  magnetic  local  time  at  t 

Word  8  =  data  error  quality  values  at  t 

Word  4-8  have  not  been  implemented  in  the  RBTM  and  are  set  to  -999.  Data  are  expected  in 
time  increasing  order. 

The  Kp  values  should  be  input  such  that  the  data  becomes  a  step  function.  For  example  on  day  2 
of  1993.0,  if  Kp  =1  for  hours  0-3  and  Kp  =  2  for  hours  3-6  the  file  should  look  like 

1.000  1993.000  2.0000  -999  -999  -999  -999  -999 

1.000  1993.000  2.1240  -999  -999  -999  -999  -999 

2.000  1993.000  2.1250  -999  -999  -999  -999  -999 

2.000  1993.000  2.2490  -999  -999  -999  -999  -999 

The  Dst  values  should  be  shifted  to  the  middle  of  the  hour.  For  example  if  Dst  is  -20.0  for 
hours  0-1  and  Dst  is  -40.0  for  hours  1-2  the  file  should  look  like 

-20.000  1993.000  2.020833  -999  -999  -999  -999  -999 
-40.000  1993.000  2.062500  -999  -999  -999  -999  -999 
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Static  Files  Required  for  Runs 


File 

Description 

boxxxxxxx.dat 

The  magnetic  field  matrices  needed  for  the 
model  run 

COORD 

Values  used  to  set  up  the  coordinate 
system 

DKTABLE 

Loss  lifetimes  for  computing  ion  loss  by 
charge  exchange 

EFCOEF 

Coefficients  fi-om  Heppner-Maynard 
model  which  are  input  to  the  electric  field 
model 

HARDY 

Coefficients  used  by  the  Hardy  electron 
precipitation  model 

lONENG 

Coefficients  used  as  input  for  the  ion 
precipitation  model 

lONNUM 

Coefficients  used  as  input  for  the  ion 
precipitation  model 
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Data  Files  Required  for  Runs 


File 

Description 

Format 

EBINS 

Number  and  value  ranges  of 
precipitating  particle  fluxes 
to  be  computed 

EBINS  (see  attached 
documentation) 

ENCHAN 

Input  values  for  the  nmnber 
of  energy  channels,  particle 
type  and  energy  channels 
for  the  model  to  simulate 

ENCHAN  (see  attached 
documentation) 

IMFBX 

Hourly  averaged  Bx 
component  of  the 
interplanetary  magnetic 
field  in  GSM  coordinates 
for  use  in  the  Kp  neural 
network 

INDATA 

IMFBY 

Hourly  averaged  By 
component  of  the 
interplanetary  magnetic 
field  in  GSM  coordinates 
for  use  in  the  Kp  neural 
network 

INDATA 

IMFBZ 

Hourly  averaged  Bz 
component  of  the 
interplanetary  magnetic 
field  in  GSM  coordinates 
for  use  in  the  Kp  neural 
network 

INDATA 

SWVEL 

Hourly  averaged  solar  wind 
velocity  values  for  use  in 
the  Kp  neural  network 

INDATA 
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RBTMIN 

Specify: 

RBTMIN  (see  attached 

1.  Start  and  stop  times 

2.  Logical  record  number 
for  restart 

3.  80  character  run 
identification 

documentation) 

4.  Sunspot  number 

5.  Output  file  prefix 

6.  Print  control  variable 

0  =  no  print 

1  =  print 

7.  Cross  polar  cap 
correction  factor  (set  to 

1.34  prior  to  DBASE4) 

8.  Kpmode 

0  =  full-up  RBTM  run 
using  all  available  data 

1  =  use  Kp  to  produce 
proxy  inputs 

9.  Forecast  mode  variable 

0  =  do  not  forecast 

1  =  forecast 
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EBINS  Format 


This  file  contains  the  number  and  value  ranges  in  eV  of  precipitating  particle  fluxes  to  be 
computed  within  the  RBTM. 


Example 


4 

30., 100. 

100..  5000. 

5000.,  15000. 
moo. ,9999999. 


Number  of  precipitating  particle  bins  (Maximum  4) 
30-  100  eV  energy  range 
100  eV  -  5  keV  energy  range 
5  keV  - 15  keV  energy  range 
>15  keV  energy 
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ENCHANFomat 


Input  file  ENCHAN  controls  how  many  electron  energy  channels  are  calculated. 
The  particle  type  is: 

1  for  electrons 


The  energy  channels  are  in  eV  with  the  energies  in  ascending  order.  If  the  energy  is  less  than 
0.65  MeV  or  greater  and  5.75  MeV  the  energy  is  reset. 

Example 

3  Number  of  energy  channels  (Maximum  37) 

1  700000.  Electrons  with  energy  0.7  MeV 
1  1000000.  Electrons  with  energy  10  MeV 
1  1500000.  Electrons  with  energy  15  MeV 


The  particle  types  are  integer;  the  particle  energies  are  floating  point. 


If  the  niimber  of  energy  channels  is  set  to  0,  the  following  default  energies  are  used: 
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1  650000. 
1  950000. 
1  1600000. 
1  2000000. 
1  2350000. 
1  2750000. 
1  3150000. 
1  3750000. 
1  4550000. 
1  5750000. 


Electrons  with  energy  0.65  MeV 
Electrons  with  energy  0.95  MeV 
Electrons  with  energy  1.6  Mev 
Electrons  with  energy  2.0  MeV 
Electrons  with  energy  2.35  MeV 
Electrons  with  energy  2.75  MeV 
Electrons  with  energy  3.15  MeV 
Electrons  with  energy  3.75  MeV 
Electrons  with  energy  4.55  MeV 
Electrons  with  energy  5.75  MeV 


If  the  number  of  energy  channels  is  set  to  0,  care  should  be  taken  with  applications  which  rely  on 
the  ENCHAN  channel  to  give  the  energy  channels  for  the  run. 


17 


RBTMIN  Format 


The  RBTMIN  file  contains  the  input  and  output  times  and  other  control  variables  for  a  particular 
RBTMrun 

Example 

Start  year,  start  day,  start  time  in  seconds 
End  year,  end  day,  end  time  in  seconds 
Beginning  record  number  (0  signifies  a  new  start) 

Up  to  80  character  run  identification 
Sunspot  number  (see  below  for  more  information) 

Output  file  prefix 
Diagnostic  print  control  variable 
0  =  no  print 
1=  print 

Cross  polar  cap  correction  factor  (set  to  1.34  for  SSIES  data  generated 
prior  to  DBASE4  which  went  online  in  July  1993) 

Kpmode 

0  =  full-up  RBTM  run  using  all  available  data 
1  =  use  Kp  to  produce  proxy  inputs 
Forecast  mode  variable  (see  below  for  more  information) 

0  =  do  not  forecast 
1  =  forecast 

The  sunspot  number  is  the  standard  published  simspot  number  and  is  available  from  the  National 
Oceanic  and  Atmospheric  Administration's  Solar  Geophysical  Data  available  fi'om  the  National 
Geophysical  Data  Center  in  Boulder,  Colorado  (see  example  attached).  It  is  also  available  at  the 
Web  page  of  the  Space  Environment  Center  in  Boulder.  For  an  event,  the  sunspot  number  can 
be  averaged.  If  the  simspot  number  is  unavailable,  the  value  of  50.0  should  be  used  as  a  default. 

To  forecast  Dst  1  hour  ahead  the  following  parameters  are  required: 

IMF  Bx,  By,  Bz  from  -3  hours  to  current  time 
Solar  Wind  Density  current  time 
Solar  Wind  Velocity  current  time 
Dst  from  -2  hours  to  current  time 

To  forecast  the  cross  polar  cap  potential  30  minutes  ahead  the  following  parameters  are  required: 

IMF  Bx,By,Bz  from  -2  hours  to  current  time 

Solar  Wind  Velocity  from  -2  hours  to  current  time 

To  forecast  the  equatorward  edge  of  the  auroral  oval  1  hour  ahead  the  following  parameters  are 
required: 

IMF  Bx,By,Bz  from  -2  hours  to  current  time 


1992  57  46800 
1992  57  48600 
0 

'30  minute  test  run' 

50.0 

'a' 

1 


1.00 

0 


0 
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Solar  Wind  Velocity  or  cross-polar  cap  potential  from  -2  hours  to  current  time 


Data  Files  Optional  for  Run 


File 

Units 

Description 

Format 

DST 

Nanotesla 

Dst  values  for  the  event 

INDATA 

EQEDGE 

Degrees 

Low  latitude  boundary  of 
the  auroral  oval  projected  to 
midnight 

INDATA 

PCP 

Kilovolts 

Polar  cap  total  potential 
drop  for  event  used  as  input 
to  the  electric  field  model 

INDATA 

SUMKP 

None 

The  sum  of  the  3-hour  Kp 
index  for  each  of  the  10 
days  preceding  the  current 
day.  These  data  are  used  by 
the  high-energy  electron 
subroutine. 

INDATA 

SWDEN 

cm"^ 

Solar  wind  density  values 
used  to  calculate  the 
standoff  distance  for  the 
event 

INDATA 

XPATT 

None  (see 

attached 

documentation) 

Polar  cap  potential  pattern 
used  as  input  to  electric 
field  model  for  the  event 

INDATA 
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Output  Files 
Standard  Output 


File 

Units 

Description 

aloct 

Radians  eastward 
from  noon 

Grid  local  time 
array 

augpar 

Various 

Augmented  data 
array  for  input 
values 

bmin 

nT 

Equatorial  magnetic 
field  strength 

bndloc 

None 

Location  of  outer 
boundary  of  detailed 
particle  traces 

colat 

Radians 

Grid  colatitude  array 

eavg 

eV 

Average  energy  in 
precipitating  energy 
bins  for  electrons 

flux 

#/cm2-s-keV 

fluxke 

eV 

flxbnd 

#/cm2-s-eV-sr 

flxsum 

ergs/cm2-s 

Binned  precipitating 
energy  flux  for 
electrons 

infofile 

None 

Information  output 
from  model  run 

ipiflx 

eV 

Average  energy  in 
precipitating  energy 
bins  for  ions 

ipieng 

ergs/cm^-s 

Binned  precipitating 
energy  flux  for  ions 

mode 

None 

Integer  array 
containing 
information  on  the 
sources  of  the  input 
data  used 

Format _ 

2-d  OUTPUT 
format _ 

1- d  OUTPUT 
format 

2- d  OUTPUT 

format _ 

1- d  OUTPUT 
format 

2- d  OUTPUT 

format _ 

5-d  OUTPUT 
format 

3- d  OUTPUT 
format 


3-d  OUTPUT 
format 


2-d  OUTPUT 
format 


5-d  OUTPUT 
format 

Ascii  text 

5-d  OUTPUT 
format 

5-d  OUTPUT 

format _ 

1-d  OUTPUT 
format 
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V 

Volts 

Electric  potential 
distribution  on  grid 
(average  of  vnrth 
and  vsoth) 

2-d  OUTPUT 
format 

vm 

(Re/nT)-2/3 

(Flux  tube 

volume)"2/3 

2-d  OUTPUT 
format 
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File 

Units 

Format 

vnrth 

Volts 

Northern 

hemisphere  electric 

potential 

distribution 

2-d  OUTPUT 
format 

vsoth 

Volts 

Southern 

hemisphere  electric 

potential 

distribution 

2-d  OUTPUT 
format 

xmin 

Re 

GSM  X  location  of 
where  field  line 
point  through  grid 
point  crosses  the 
equatorial  (B-field 
minimum)  plane 

2-d  OUTPUT 
format 

ymin 

Re 

GSM  Y  location  of 
where  field  line 
point  through  grid 
point  crosses  the 
equatorial  (B-field 
minimum)  plane 

2-d  OUTPUT 
format 

zmin 

Re 

GSM  Z  location  of 
where  field  line 
going  through  grid 
point  crosses  the 
equatorial  (B-field 
minimum)  plane 

2-d  OUTPUT 
format 
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Standard  OUTPUT  Interface 


The  standard  call  to  subroutine  OUTPUT  is: 

CALL  OUTPUT(LUN,  IRECMX,  ID,  RID,  CHID,  ARRAY,  IDIM,  JDM,  KDM,  IMAX, 
JMAX,  KMAX,  PREFIX,  FILNAM) 

LUN  Logical  Unit  number  on  which  to  write  file 

IRECMX  Record  Number  to  begin  write  function 

ID  Integer  header  vector 

ID(1)  Year 

ID(2)  Day 

ID(3)  Hour 

ID(4)  Minutes 

ID(5)  Seconds 

ID(6)  Presently  unused 

ID(7)  Presently  imused 

ID(8)  First  dimension  of  output  array 

ID(9)  Second  dimension  of  output  array 

ID(  1 0)  Third  dimension  of  output  array 

ID(ll)  Time  index  L 

ID(1 2-20)  Presently  unused 

RID  Real  header  vector 

RID(l)  Time  tag  (seconds) 

RID(2)  Kp 

RID(3)  Polar  cap  potential  drop  (kV) 

RID(4)  Time  derivative  of  location  of  low-latitude 
edge  of  auroral  oval  (degrees/hour) 

RID(5-20)  Presently  unused 

CHID  Character  header  string  (up  to  80  characters  long)  read  in  from  RBTMIN. 

ARRAY  Array  to  be  written 

EDIM  I  dimension  of  array 

JDIM  J  dimension  of  array 

KDIM  K  dimension  of  array 

IMAX  Actual  I  maximum  in  output  ARRAY 

IMAX  Actual  J  maximum  in  output  ARRAY 

KMAX  Actual  K  maximum  in  output  ARRAY 

PREFIX  Run  identification  character  from  RBTMIN 

FILNAM  File  to  be  written 
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1-d  OUTPUT  Format 


These  files  are  written  by  SUBROUTINE  OUTPUT  and  have  the  same  structure: 

Record  IID,  RID,  CHID,  ARRAY  for  timestep  1 
Record  2ID,  RID,  CHID,  ARRAY  for  timestep  2 


Record  nID,  RID,  CHID,  ARRAY  for  timestep  n 

For  the  following  1-d  OUTPUT  files  the  only  parameter  that  changes  is  the  record  length 
because  of  the  different  ARRAY  sizes: 

augpar  Dimension  NAUGEL  (currently  set  to  28) 

The  augpar  array  contains: 


augpar(l) 

Year 

augpar(2) 

Day 

augpar(3) 

Seconds  of  day 

augpar(4) 

Minute 

augpar(5) 

Seconds 

augpar(6) 

Kp 

augpar(7) 

Dst  (nT) 

augpar(8) 

Equatorward  edge  of  the  auroral  oval  at  midmght  (degrees) 

augpar(9-ll) 

Presently  unused 

augpar(12) 

Bx  component  of  the  interplanetary  magnetic  field  (nT) 

augpar(13) 

By  component  of  the  interplanetary  magnetic  field  (nT) 

augpar(14) 

Bz  component  of  the  interplanetary  magnetic  field  (nT) 

augpar(15) 

Magnetotail  collapse  parameter 

augpar(16-20) 

Presently  unused 

augpar(21) 

Solar  wind  velocity  (km/s) 

augpar(22) 

Solar  wind  density  (cm"3) 

augpar(23) 

Presently  unused 

augpar(24) 

Standoff  distance 

augpar(25) 

Cross  polar  cap  potential  drop  (kV) 

augpar(26) 

Electric  field  pattern  type 

augpar(27) 

Time  rate  of  change  of  Dst  (nT/hour) 

augpar(28) 

Time  rate  of  change  of  the  equatorward  edge  of  the  auroral  oval  at 
midnight  (degrees/hour) 

bndloc  Dimension  JDIM  (currently  set  to  5 1) 

The  bndloc  array  contains  the  floating  point  i  value  of  the  outer  boxmdary  of  the  particle 
traces  for  each  local  time  0)  gridpoint. 
mode  Dimension  NAUGEL  (currently  set  to  28) 

For  each  of  the  28  augpar  variables  given  above,  the  mode  array  contains  information  on  where 
the  input  data  was  obtain  firom: 
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0  Default  (Kp-based)  front-end  model  used 

1  Data  from  input  data  stream 

2  Forecast  module  value  used 

3  Persistence  value  used 

4  Interpolated  value  between  forecast  and  data  stream  value 

5  Interpolated  value  between  data  stream  and  default  values 

6  Interpolated  between  forecast  and  default  values 
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2-d  OUTPUT  Format 


These  files  are  written  by  SUBROUTINE  OUTPUT  and  have  the  same  structure  as  the  1-d 
OUTPUT  format: 

Record  IID,  RID,  CHID,  ARRAY  for  timestep  1 
Record  2ID,  RID,  CHID,  ARRAY  for  timestep  2 

Record  nID,  RID,  CHID,  ARRAY  for  timestep  n 
The  ID,  RID  and  CHID  also  have  the  same  structure. 

The  following  arrays  are  dimensioned  IDIM  (currently  set  to  62)  by  JDIM  (currently  set  to  62) 
aloct 
bmin 
colat 

V 

vm 

vnrth 

vsoth 

xmin 

ymin 

zmin 

The  following  array  is  dimensioned  JDIM  (currently  set  to  51)  by  KDIM  (currently  set  to  37) 
flxbnd 
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3-d  OUTPUT  format 


These  files  are  written  by  SUBROUTINE  OUTPUT  and  are  used  by  the  flux  and  fluxke  output 
files.  For  k  energy  channels  the  format  is: 

Record  IID,  RID,  CHID,  ARRAY  for  energy  1  for  timestep  1 
Record  2ID,  RID,  CHID,  ARRAY  for  energy  2  for  timestep  1 

Record  kID,  RID,  CHID,  ARRAY  for  energy  k  for  timestep  1 
Record  k+1  ID,  RID,  CHID,  ARRAY  for  energy  1  for  timestep  2 
Record  k+2  ID,  RID,  CHID,  ARRAY  for  energy  2  for  timestep  2 

Record  k+n  ID,  RID,  CHID,  ARRAY  for  energy  k  for  timestep  n 
The  ID,  RID  and  CHID  also  have  the  same  structure. 

The  flux  arrays  are  dimensioned  IDIM  (currently  set  to  62)  by  JDIM  (currently  set  to  62)  by 
KDIM  (maximum  37) 


5-d  OUTPUT  Format 


These  files  are  written  by  SUBROUTINE  OUTPUT  and  are  used  by  the  precipitating  number 
and  average  energy  electron  and  ion  output  files.  The  number  of  bins  is  set  in  file  EBINS,  with 
the  maximum  number  of  4 

Record  1  ID,  RID,  CHID,  ARRAY  for  bin  1  for  timestep  1 
Record  2ID,  RID,  CHID,  ARRAY  for  bin  2  for  timestep  1 
Record  3ID,  RID,  CHID,  ARRAY  for  bin  3  for  timestep  1 
Record  4ID,  RID,  CHID,  ARRAY  for  bin  4  for  timestep  1 
Record  SID,  RID,  CHID,  ARRAY  for  bin  1  for  timestep  2 


Record  n+4  ID,  RID,  CHID,  ARRAY  for  bin  4  for  timestep  n 
The  ID,  RED  and  CHID  also  have  the  same  structure. 

The  arrays  are  dimensioned  IDIM  (currently  set  to  62)  by  JDIM  (currently  set  to  62)  by 
KBNDM  (currently  set  to  4)  by  ITMDIM  (currently  set  to  50) 

The  following  files  use  this  format: 
flxsum 
eavg 
ipiflx 
ipieng 
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APCAL 


a.  Function  -  Function  to  calculate  Apl  5  from  neural  net  coefficients. 

b.  Input - 

FKPIN  Vector  of  Kp  values  for  current  time  step  to  timestep  15  (3.5  hours) 

c.  Processing  -  Uses  neural  network  coefficients  to  calculate  Apl 5. 

d.  Output - 

APCAL  Apl  5  value  for  current  time  step 
APNN 

a.  Function  -  Subroutine  to  calculate  Apl 5  from  neural  net  coefficients. 

b.  Input - 

ITMMAX  Maximum  number  of  time  labels 
NELTS  Number  of  elements  in  output  data  array 
ITMDIM  Maximum  number  of  time  steps  per  run 
NAUGEL  Number  of  elements  in  the  augmented  input  array 
PREFIX  Prefix  for  current  run 

c.  Processing  -  following  are  the  subroutines  called  and  their  major  fimctions 

APCAL  Uses  neural  network  coefficients  to  calculate  Apl 5. 

This  routine  reads  in  the  necessary  input  data  for  the  Apl 5  neural  network  and  calls  the  fimction 
APCAL  to  calculate  the  Apl 5  value  for  each  timestep  in  the  run.  The  neural  network  requires 
the  current  Kp  and  the  Kp  value  for  15  previous  time  steps.  If  this  is  a  new  run,  all  values  are  set 
using  the  current  Kp  value.  The  calculated  Apl 5  value  is  put  into  PARRAY(20,itm). 

d.  Output - 

PARRAY  Data  array  of  interpolated  data  values 
MODE  Integer  variable  denoting  input  data  sources 
0  Default  (Kp-based)  front-end  model  used 

1  Data  from  input  data  stream 

2  Forecast  module  value  used 

3  Persistence  value  used 

4  Interpolated  value  between  forecast  and  data  stream  value 

5  Interpolated  value  between  data  stream  and  default  values 

6  Interpolated  between  forecast  and  default  values 


29 


BNDCRS 


a.  Function  -  Set  up  boundary  flux  at  2.5  and  6.75  Re  from  CRRESELE  model.  Set  up 

geosynchronous  orbit  flux  at  6.6  Re. 

b.  Input - 

LATDIM  Number  of  latitudinal  grid  lines 

LTDIM  Number  of  local  time  grid  lines 

ITMDIM  Maximum  number  of  times  steps  per  run 

ITMMAX  Maximum  number  of  time  labels 

lEDIM  Maximum  number  of  energy  chaimels  per  run 

lEMAX  Number  of  energy  chaimels  in  current  run 

NAUGEL  Number  of  elements  in  the  augmented  input  array 

FLXCRS  Flux  by  energy  and  L  shell  from  the  CRRESELE  model  (electrons/cm^-s-keV) 
XKE  Energy  channels  used  in  current  run  (MeV) 

AUGPAR  Augmented  array  of  input  values 
R  RBTM  grid  in  Re  for  current  time  step 

COLAT  Grid  colatitude  array  (radians) 

LL  Number  of  current  time  step 

c.  Processing  -  the  following  are  the  subroutines  called  and  their  major  functions: 

FNDBI  Determines  BI  value  corresponding  to  a  given  R  for  a  specified  J  line 

GEOFLX  Subroutine  to  read  file  of  geosynchronous  fluxes  and  place  in  RBTM  grid. 

Reads  file  of  geosynchronous  flux  and  locates  them  in  the  RBTM  grid. 

d.  Output  - 

FLXCRB  Flux  at  2.5  Re  in  RBTM  grid  for  all  energies  (electrons/cm^-s-keV) 

BNCREQ  I  value  of  point  at  2.5  Re 

FLXPL  Flux  at  6.75  Re  in  RBTM  grid  for  all  energies  (electrons/cm2-s-keV) 

BNDPL  I  value  of  point  at  6.75  Re 

FLXGEO  Flux  at  geosynchronous  orbit  for  all  local  time  for  current  time  step  for  given 
energy  (electrons/cm^-s-keV) 

GEOFLX 

a.  Function  -  Subroutine  to  read  file  of  geosynchronous  fluxes  and  place  in  RBTM  grid. 

b.  Input  - 

LATDIM  Number  of  latitudinal  grid  lines 

LTDIM  Number  of  local  time  grid  lines 

ITMDIM  Maximum  number  of  time  steps  per  run 
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ffiDM 

NAUGEL 

ITM 

lEMAX 

AUGPAR 


Maximum  number  of  energy  channels  per  run 
Number  of  elements  in  the  augmented  input  array 
Current  time  step 

Number  of  energy  channels  in  current  run 
Augmented  data  array  for  input  values 


c.  Processing - 

d.  Output - 

FLXGEO  Flux  at  geosynchronous  orbit  for  all  local  time  for  current  time  step 
(electrons/cm^-s-keV) 


INTCRS 


a.  Function  -  Set  up  initial  flux  distribution  using  CRRESELE  model. 


b.  Input - 


LATDIM 

LTDIM 

lEDIM 

lEMAX 

ITMDIM 

R 

XKE 

XMOD 


Number  of  latitudinal  grid  lines 
Number  of  local  time  grid  lines 
Maximxun  number  of  energy  channels  per  run 
Number  of  energy  channels  in  current  run 
Maximum  number  of  times  steps  per  run 
RBTM  grid  in  Re  for  current  time  step 
Energy  chaimels  used  in  current  run  (MeV) 
Apl5  value  for  initial  time  step 


c.  Processing  -  Reads  file  of  flux  and  L-shell  data  from  CRESSELE  model  data  files  and  places 
them  onto  RBTM  grid. 


d.  Output - 

FLXCRS  Flux  values  from  CRRESELE  model  by  energy  by  L-shell 
(electrons/cm^-s-keV) 

BIN  L  values  used  in  CRRESELE  model. 

FLXBEG  Initial  flux  values  for  each  energy  on  RBTM  grid  (electrons/cm^-s-keV) 


KPNN 

a.  Frmction  -  Calculate  Kp  using  neural  network  coefficients  and  enter  into  DARRY. 

b.  Input  - 

NDIM  Dimension  of  DARRY  giving  maximum  number  of  data  records 

DARRY  Observational  data  array 

NUMNUM  Number  of  observational  data  elements 
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ISTART  Rim  start  time  (year, day, seconds) 

c.  Processing  -  The  following  are  the  subroutines  called  and  their  major  functions 

FORTIM  Performs  model  time  conversion 
KPNNCON  Kp  calculation  controller 

If  this  is  the  1st  time  step,  the  forecast  Kp  value  is  put  in  the  T+1  hour  position  in  DARRY  and 
the  T,  T+15  minutes,  T+30  minutes  and  T+45  minutes  values  are  backfilled  with  the  same  Kp. 

d.  Output - 

DARRY  Observational  data  array  with  neural  net  Kp  entered 
KPNNCON 

a.  Function  -  Kp  neural  net  controller 

b.  Input - 

YEAR  Year 

FDAY  Decimal  day 

c.  Processing  -  Following  are  the  subroutines  called  and  their  major  function: 

GETHIS  Parameter  history  acquisition  routine  for  the  Kp  neural  net 
KPNNMOD  Kp  neural  network  model 

d.  Output  - 

Kp  Neural  net  calculated  Kp 

lERR  Error  return 

=0  Forecast  complete 

=1  No  forecast  because  of  missing  data 


KPNNMOD 

a.  Function  -  Function  to  calculate  Kp  from  neural  net  coefficients. 

b.  Input - 

B  IMF  hourly  average  B-field  magnitude  for  time  T  and  time  T-1  hour  (nT) 

By  IMF  hourly  average  of  the  Y-component  for  time  T  and  time  T-1  hour  (nT) 

Bz  IMF  hourly  average  of  the  Z-component  for  time  T  and  time  T-1  hour  (nT) 

SWVEL  Solar  wind  speed  homly  average  for  time  T  and  T-1  hour  (Km/s) 
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c.  Processing  -  Uses  neural  network  coefficients  to  calculate  Kp. 

d.  Output - 

FORVAL  Kp  value  for  hour  T+1 
lERR  Error  flag 

=1  Problem  in  data  normalization  (value  is  >  1  or  <  0)  the  value  is 
reset  (to  0  or  1)  and  the  neural  network  continues 
=2  Calculate  Kp  is  <  0  or  >  9 


PARGEN 

a.  Function  -  Subroutine  to  obtain  values  from  the  environmental  data  base  and  interpolate  or 
extrapolate  as  appropriate  to  provide  data  at  a  normalized  time. 

b.  Input - 


ISTART 

lESfC 

lEND 

ITMDIM 

NELTS 

ITMMAX 

TIMTAG 

NAUGEL 

XIP 

XPCP 

NIP 

NPCP 

PREFIX 


Start  time  (year, day, seconds) 

Increment  time  (year, day, seconds) 

Stop  time  (year,day,seconds) 

Maximum  number  of  times  steps  per  run 
Number  of  elements  in  output  data  array 
Maximum  number  of  time  labels 

Vector  giving  times  at  which  E  and  B  parameters  are  calculated 

Number  of  elements  in  the  augmented  input  array 

Saved  IPATT  DARRY 

Saved  PCP  DARRY 

Number  of  data  elements  in  XIP 

Niunber  of  data  elements  in  XPCP 

Single  character  prefix  for  current  run 


c.  Processing  -  The  following  are  the  subroutines  called  and  their  major  functions: 


APNN 

INDATA 

DTACHK 

SMOOTH 

DTXIPT 

DTNTRP 

TIMINC 

TCONV3 


Subroutine  to  calculate  Apl5  from  neural  net  coefficients 
Subroutine  to  read  in  data  for  the  MSFM 
Subroutine  to  check  for  input  data  values  out  of  range 
Subroutine  to  extract  data  from  input  array 
Subroutine  to  return  interpolated  polar  cap  patterns 
Subroutine  to  return  interpolated  data  values 

Subroutine  to  increment  time  for  next  electric  and  magnetic  field  record 
Function  subprogram  to  convert  to  standard  program  representation  of  time. 


This  subroutine  retrieves  observational  data  from  the  environmental  data  base  and  interpolates 
as  necessary  to  have  input  data  form  the  entire  run. 
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d.  Output  - 


PARRAY  Data  array  of  interpolated  data  values 
MODE  Integer  variable  denoting  input  data  sources 
0  Default  (Kp-based)  ffont-end  model  used 

1  Data  from  input  data  stream 

2  Forecast  module  value  used 

3  Persistence  value  used 

4  Interpolated  value  between  forecast  and  data  stream  value 

5  Interpolated  value  between  data  stream  and  default  values 

6  Interpolated  between  forecast  and  default  values 

RADBLT 

a.  Function  -  Calculate  Radiation  Belt  electron  flux  using  fully  adiabatic  method. 


b.  Input - 


lEMAX 

Cl 

FLUXl 

XKE 

FNDPLl 

BNDPLl 

FBNDEQl 

BNDEQl 

B1 

B2 

C2 

FNDPL2 

BNDPL2 

FBNDEQ2 

BNDEQ2 

R 

FLXGEO 


Niunber  of  energy  channels  in  current  run 

Grid  colatitude  array  for  previous  time  step  (radians) 

Flux  by  energy  on  RBTM  grid  for  previous  time  step  (electrons/cm^-s/keV) 
Energy  channels  used  in  current  run  (MeV) 

Flux  at  6.75  Re  in  RBTM  grid  for  all  energies  for  previous  time  step 
(electrons/cm^-s-keV) 

I  value  of  point  at  6.75  Re  for  previous  time  step 

Flux  at  2.5  Re  in  RBTM  grid  for  all  energies  for  previous  time  step 

(electrons/cm^-s-keV) 

I  value  of  point  at  2.5  Re  for  previous  time  step 

Magnetic  field  magnitude  at  equatorial  plane  for  previous  time  step  (nT) 
Magnetic  field  magnitude  at  equatorial  plane  for  current  time  step  (nT) 

Grid  colatitude  array  for  current  time  step  (radians) 

Flux  at  6.75  Re  in  RBTM  grid  for  all  energies  for  current  time  step 

(electrons/cm^-s-keV) 

I  value  of  point  at  6.75  Re  for  current  time  step 

Flux  at  2.5  Re  in  RBTM  grid  for  all  energies  for  current  time  step 

(electrons/cm^-s-keV) 

I  value  of  point  at  2.5  Re  for  current  time  step 
RBTM  grid  in  Re  for  current  time  step 

Flux  at  geosynchronous  orbit  for  all  local  time  for  current  time  step 
for  given  energy  (electrons/cm^-s-keV) 


c.  Processing  -  The  following  are  the  subroutinies  called  and  their  major  functions: 

FLXCALC  Calculates  energy  and  flux  change  on  each  L  shell  and  places  them 
onto  RBTM  grid 


> 


FRL  Function  to  calculate  phi  (3rd  adiabatic  invariant) 

G3NTRP  Generic  3d  interpolation  function 

SETUP  Calculates  L  shells  and  location  from  magnetic  field  of  previous  time  step 

SHLDFT  Calculates  L  shells  and  location  from  magnetic  field  of  current  time  step 
YFIT  Function  to  do  linear  interpolation 

d.  Output  - 

FLUX2  Flux  by  energy  on  RBTM  grid  for  current  time  step  (electrons/cm^-s/keV) 
BNDLOC  Outer  boundary  set  at  8.0  Re 
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